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ABSTRACT 

A finite number of harmonic oscillators coupled to infinitely many environment 
oscillators is fundamental to the problem of understanding quantum dissipation of 
a small system immersed in a large environment. Exact operator solution as a 
function of time is given to this problem, by using diagonalized dynamical variable 
of the entire system, the small system plus the environment. The decay law of 
prepared initial configuration is worked out in greatest detail. A clear separation of 
the exponential- and the power-law decay period is made possible by our method. 
Behavior of physical quantities at asymptotically late times can be understood in 
terms of the overlap probability of the system variable with the diagonal variable of 
the entire system. 
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How a small system behaves under influence of a larger environment has been 
investigated in a variety of approaches. The simplest, yet the most fundamental 
model of quantum dissipation is harmonic oscillator coupled to infinitely many os- 
cillators that make up a bath in a mixed state. Two powerful methods to analyze 
this problem are the quantum Langevin equation [||], and the path integral approach 
im]. In our previous work p| we used the path integral method for the linear open 
system, and derived the reduced density matrix of the system, by integrating out 
environment variables. In the present work we employ a direct operator approach. 
Our new approach has an advantage, giving explicit operator solution and thereby 
extends previous results to arbitrary initial states including correlation between 
the system and the environment, and to non-thermal environments. 

Consider a harmonic oscillator denoted by its coordinate variable g, interacting 
with bath oscillators Q{uj)^s, labeled by the continuous frequency u. The Hamil- 
tonian of the entire system is assumed to be given by a sum of the kinetic term, 
I (p^ + / dujP'^{ijj) ) and the potential term, 

written in the matrix notation. c{uj) describes the coupling of the variable g, which 
we regard as the small system in question, to the environment variable Q{uj). The 
potential part may be decomposed into the two, V = Vo + V , where Vo consists of 
the environment part alone, 

n \ 

(2) 

{S,\S,)V,, V^ = V,, (3) 

(4) 

duj2ur{u!), (5) 

where Uc is the smallest value of the environment spectrum. r{uj) = c^(u;)/(2co') 
characterizes the response of environment on the system. Our results in the present 
work hold for any r{uj) obeying certain convergence property. 

A typical question one frequently asks with regard to the system behavior is how 
a pure state of the q system evolves in time in environment described by a mixed 
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state such as the thermal one of Q{uj)^s. We shall answer this by giving explicit 
time dependent operator solution for q{t) ,p{t) that is written in terms of initial 
values, qi ,Pi , Qiiyj) , Piiuj) . With this one can clearly express correlators such as 
{q{ti)q{t2)) in terms of the distribution of these initial values in any state. Closest to 
our present work is the classical work of Ullersma [0], with an important difference of 
the presence of a gap in the environment spectrum, which yields different behaviors 
of the correlators. The gapless case is applied for instance to phonons in medium, 
but in the case of our interest such as the unstable particle decay in medium, a gap 
exists if the mass of the daughter particle is finite. 

The model itself, a harmonic system coupled to an infinite number of harmonic 
oscillators, frequently appears in many idealized physical problems. We utilize in the 
present investigation exact solution of the Hamiltonian eigenvalue problem to this 
system, which might also be useful in other contexts. Essentially the same solution of 
diagonalization is given in ref 0] in a related model, but we use a modified form with 
an extra phase factor attached, which makes our subsequent derivation transparant. 

We first solve a related problem by seeking the eigenvector |\l'(co')) diagonalizing 
the matrix V, using the unperturbed base \uj), ( Vq — cu^ ) \uj) = 0; 

We have introduced the out-state boundary condition: The boundary condition is 
determined by the prescription of the limit, ilim^_,o+ e = iO~^ , which extracts 
well-behaved solutions at t — > +oo, as in the Lippman-Schwinger equation. This 
choice of the boundary condition is superficial and not dictated, to the extent that 
our problem is diagonalization by orthogonal transformation of hermitian oscillator 
variables. But we shall soon see how the solution with this form is useful to solving 
our problem. 

Our ansatz of the eigenvalue problem (P) assumes that there exists no isolated 
spectrum outside the continuos one. An isolated spectrum, in the context of the 
problem raised here, corresponds to a bound, stable state, whose presence would 
imply that an initially prepared configuration of the system does not completely 
decay by interaction with environment. We exclude this case for our consideration. 
Let us clarify under what condition the spectrum without any isolated spectrum is 
realized. The equation that governs the real part of the isolated energy eigenvalue 
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squared A is given by 

/(A) = A, fiX)^ujl + V r dJ-^^. (7) 

From the monotonic behavior of /(A) in the X < uj^ region one readily derives the 
condition of no zero in this region: 

fi^c') > , (8) 

with ~ indicating the hmit from below. Consistency further requires uoq > ujc- When 
this condition (^) is obeyed, the system energy ojq moves, with increase of the cou- 
pling c(co'), into the second Riemann sheet of the complex uj plane having a cut at 
uj > uJc- The imaginary part of this location gives the decay rate of any initial 
configuration of the system, as will be shown later. 

From the ansatz of above we get equations for q = {Si\^{ijj)) , (z = 1 , 2) 

and thus eigenvectors. 




(9) 



We introduced the analytic function F{z) 0; 

1 duo riuo) , . 

with extended r{—uj) = — r(u;) for u < 0. This function has cuts along the real 
axis, u > Uc and u < —uJc- We assume Uc > such that there is a gap between 
the two cuts. The following discontinuity relation holds; 

F{uj + iO+) - F{uj - iO+) = i2nr{u)F{uj + iO+)F(cu - iO+) = i27iH{u) , (11) 

along this cut. 

The previous function /(A) is related to this function by 

f{X)-X = ^F-\VX + iO+). (12) 

An important property of F{z) is that in the cut 2;— plane this function is regular 
except on the real axis, and that it can be analytically continued into the other 
Riemann sheets by the discontinuity formula. In the second sheet F{z) can be 
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shown to have simple poles corresponding to unstable states. The location of these 
poles is given by 

-ujI-27iG{z) + 2mr{z) = 0. (13) 

Using the discontinuity formula, we may determine the the normalization, leading 
to |A^(co')p = 1 . Hence the overlap probability of the two vectors is 

\{0\-q/{uj))\'^ = 2ujH{uj) , (14) 

r(co') 

^^""^ " (a;2-a;2-n(a;))2 + (7rr(a;))2 ' ^^^^ 

with 

r°° , r{uj') 



n(cj) = V / dJ ^ ' . (16) 
J-oo CO — to' 

The eigenvector is then given by 
l^icu)) = \uj) - c{uj)F{uj - zO+) |0) + c{uj)F{uj - ^0+) / dJ—-^-^^^—— \uj') . (17) 



u;2 + iO+ 



It might seem odd due to an extra system degree of freedom, but it is possible 
to get inverted transformation to express the original base in terms of the diagonal 
base thus derived; 



oo 



|0) = -/ dcJc(cj)F*(cu-zO+) l^(cu)), (1^ 



\u;) = \^{u;)) + c{u;) / du' c{u;') ^ ' ' \^{u;')). (19) 

Jujc UJ — UJ — iU 

How are these solutions for the state vectors related to our problem of the operator 
diagonalization? The crucial fact that can be verified straightforwardly is that the 
complex coefficients for the out-state solution (|l^) has a common phase of e*'^'-'^-' = 
F{uj — iO^)/\ F{uj — iO^) I . This means that after taking out the phase factor e*'^*-'^'' 
diagonalization via real orthogonal matix becomes possible, leading to hermitian 
diagonal variables Q{uj). In the following computations we find it more convenient 
to use Q{uj) = e*'^^'^) Qiuj) by retaining these u dependent phases. 

With these remarks it should be evident that the canonical transformation be- 
tween the two descriptions exists; 

Qiuj) =Q{uj)- c{uj)F{uj - zO+) (q - r duj' ''^'^'\ Q{uj') ] , (20) 



UJ 



+ iO- 

c{uj') 



P{uj) = P{u) - c{uj)F{u - ^0^) Ip-J^ duj' P{^') , (21) 
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/■oo 



poo _ 

- / dujc{uj)F*{uj -iQ+)Q{uj) , (22) 

roo 

p = - dujc{uj)F*{uj-iO+)P{uj), (23) 

/■°° c(uj') F*(uj' - iO^) ~ 
Q{u) = Q{u) + c{u:) du' "^^^^ ^ ^ g(c.O , (24) 



c{uo')F*iuj' -i^^ 



P{uj) = Pioo) + cioo) / rfu;' P{oo') . (25) 

Assuming the canonical commutation for the original variables, one can verify that 
diagonal variables obey the correct form of the commutation relation; 
[ Q{uj) , Pioj') ] = i 5{uj — oj' ) etc. It can be proved that with the specified phases of 
Q{uj) ,P{uj) the original variables are all hermitian, as required. 

Operator diagonalization thus presented extends the work of 0] only slightly: 
their model is simpler, with only one dimensional projection operator involved in 
the interaction V', while ours has two dimensional projections. It should however be 
clear that the diagonalization is possible with any finite number of projections, with 
more complexities. 

With this diagonal variable it is easy to write Heisenberg operator solutions at 
any time t in terms of the initial operator values; for instance 

/■oo _ 

q{t) = - duJc{uj)F*{uj -iO+)Q{u} ,t) , (26) 



Q{uj,t) =cos{ujt)Qi{uj) + ^^^^^^Pi{uj). (27) 

u 

The initial values Qi , Pi are then rewritten in terms of the original variables. After 
some straightforward calculation one finds that 



q{t) = P^9{t) + q^ m - duj [f\uj- 20+) e~ + (h.c.) ) 



/•oo /-oo / Jriuj') . ^ \ 

/ du' / duH{uj) / \ '.^^ e"''^'h{u') + (h.c.) , (28) 

J Uc J —OO \ ^ — ~r lO I 

roo 

g{t) = 2 du H{uj) sm{iut) , (29) 

JiUr 



where 



and bi{iu) = ( y/uj Q{uj)+iP{uj) / ^fuj ) / v2 is the annihilation operator for environment 
harmonic oscillators. 

The formula for qit) and p(t) is much simplified by introducing 

h{u),t) = f drg{T)e-''^\ (30) 
Jo 
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which is shown to be equal to 

/oo g-j(<^-a;')t 
du^' ^(^0 • (31) 
-oo iu — UJ — LU ^ 

With this function one finds that 

q{t)=p,g{t) + q,g{t)- H duo JHj^) (h* ,t) e-'^' h{uj) + ih.c.)) . (32) 
Similarly with the aid of 

kiu,t)= t dr g(T)e- ''^^ = g(t)e-'''' + iuohiuo , t) , (33) 

JO 

the momentum operator is given by 

p{t)=p,g{t) + (li 9{t)- r dujJ^) (k*{uj,t)e-'^'bi{uj) + {h.c.)) . (34) 

The basic formula (|28|) was also derived by UUersma , who however did not ex- 
ploit these solutions in generality they really deserve. (The resonance approximation 
or the local friction approximation is often made in deriving physical consequences 
in that work.) More importantly, the work of ref 0] assumed Uc = (the gapless 
case). This special case does neither produce the power-law decay nor the enhanced 
remnant as will be discussed subsequently. 

It would be useful to clarify the physical significance of the function g{t) which 
is so basic in these formulas. For this purpose let us derive a semi-classical equation 
of motion by taking the ensemble average over the initial environment that has no 
expectation value, {Qi^uj)) = {Pi{uj)) = 0. The expectation value of the system 
variable is then written in terms of initial expectation values of the system variable: 

{P{t)) = {p^) m + (g.) 9 it) , {q{t)) = {p,} g{t) + (g,) g{t) . (35) 

Thus in the semi-classical picture the function g{t) directly controls time evolution. 
Furthermore the initial values may be eliminated to give a closed equation of motion 
for the system variable: 

^ = -n(if(,>-C(«){p>, ^ = {P>, (36) 

= = (37) 

g 9 -g^ g 9 -g^ 
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The quantities in this equation give the time dependent friction (C(t)) and the time 
dependent frequency squared {Q'^(t)), incorporating environmental effects. The ran- 
dom force that appears in the quantum Langevin equation is missing here, because 
we dropped the environment variable in taking the ensemble average. Thus g{t) 
describes an average behavior of the system variable disregarding the random force 
from environment. 

Either g{t) or g{t) can be shown to follow integro-differential equation of the 
following form {y = g ot g), 



a/(r) is the imaginary part of the kernel function that appears in the influence 



grated effects of the environment, the non-local friction in general along with the 
frequency shift. The two functions, g{t) and g{t), differ in their boundary condi- 
tions: g{0) = , g{0) = 1 . The same differential equation also appeared in the path 
integral approach 0, [|]. 

The system variable has been determined in terms of the initial operator values 
of both the system and the environment variables. Dependence on the system initial 
values Pi , qi are given by the function g(t) as defined by (PP]), which decreases first 
exponentially and then finally by an inverse power of time, as can be seen in the 
following way. Using the discontinuity formula, one may rewrite the u integration 
containing H{uj) = ( F{uj + iO^) — F{uj — iO+) ) / (27ri) along the real axis into the 
F{z) integration, the complex z running both slightly above and below the cuts. 
The factor sin(a;t) is replaced by $56^*"^* in this procedure. A half of this complex 
contour can be deformed into the second sheet, and one thereby encounters simples 
poles in the second sheet, as shown in Fig.l. We assume for simplicity that there 
exists only a single pole in the nearby second sheet. The intregral for g{t) may then 
be expressed as the sum of the pole contribution (at z = zq with '^zq < 0) in the 
second sheet and the contribution parallel to the imaginary axis passing through 
z = Uc, both in the first (I) and in the second (II) sheet 0: 




(38) 



(39) 



functional 0. Its significance on the system behavior is that it incorporates inte- 




e 



TT 



/■oo 

/ dye-y^ {Fi{uj, + iy)-Fn{uj, + iy)) , (40) 
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with K 



As seen from this formula, the pole contribution given by the first term describes 
the exponential decay which usually lasts very long during the most important phase 
of the decay period, while the rest of contribution gives the power law decay at very 
late times; oc The power — a — 1 is related to the threshold behavior of the 

response weight, r{uj) oc (tu — Uc)" ■ a > — 1 is required as a consistency of this 
approach; convergence of the uj integration. An example of g{t), computed by the 
formula (ppj), is shown in Fig.2, which clearly exhibits both the exponential and the 
power law behavior. The exponential period can be fitted precisely by the first pole 
term of eg . (^OD including the phase of oscillation, except at very early times where 
quantum mechanical effect gives the non-decay probability decreasing as 1 — 0[t^] . 
The exponent factor — Qzq > gives the decay rate of any initial configuration, and 
is related to the constant friction coefficient i] in the usual treatment of dissipation, 
namely, when quantum dissipation is assumed local in time. This can be shown by 
computing the time dependent friction C{t) in eq.(0): when the non-pole term is 
neglected in eq. (^OD , 



with an approximate form of g{t) = ^ sin(a;t) e"'''* . Both at very early and very 
late times these functions deviate from the behavior in the pole approximation: for 
instance C(t) = 0[t^] as t — >• 0+. At late times both of Q'^it) and C(t) decrease with 
powers of time. 

Thus at asymptotically late times the initial value dependence disappears leading 

to 



C(t) = 27, n\t) 




(41) 



lit) ^ - f 



'OO 



duo^H^) [F*{uj-iO+)e-''^'bi{uj) + (h.c.)) 



(42) 




(^)-(h.c.)), (43) 



due to h{uj , oo) = F{uj — iO^), which may be proved from (^). The system evolution 
is thus governed by the initial distribution of the environment variables, bi{uj) , bl^u), 
with the following probability functions; 



r{uj)\h{uj,t)\^ forg(t) 



r{Lj)\k{LJ ,t)\'^ forp(t). 



(44) 



that have limits of 



r{uj)\h{uj,oo)\'^ = H{uj) 



(45) 
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riuj)\k{uj,oo)\' = uj''H{uj) = || (0|^H) |^ (46) 

Note that H{uj), hence the overlap probabihty | (0|\[^(ti;)) p has a universal character, 
determined by general properties of the system and the environment alone, irrespec- 
tive of their particular initial states. The overlap probability between the system q 
variable and the diagonal Q{uj) respects the unitarity relation, 

/■oo poo 

/ du\{0\^{Lj))\^ = duj2uH{uj) = 1, (47) 

which can be proved for any r{uj). This indicates that the final state is described in 
terms of Q ,P variables alone. 

The important function H{u!) takes almost the form of Breit-Wigner function 



having a complex pole at = z of eq.(|T3| ), especially in the weak coupling limit. 



c{uj) 0. This leads to the approximation to neglect non-pole terms in (40) for 
practical applications. We however warn that some important physics is missed in 
this approximation as will be made clear later. 

Various physical quantities of the system under question can be computed using 
our formulas written in terms of the initial values of the system and the environment. 
We consider the situation in which the environment is described by some mixed state. 
Let us assume for simplicity no correlation between the system and the environment 
at an initial time, and furthermore that the initial environment density matrix is 
diagonal with the Hamiltonian or the number operator of each u, ni{uj) = bl{uj)bi{Lj) , 
such that {hi{ijj)hi{ijj')) = etc. Various correlation functions are then evaluated; 
for instance, 



{(l{ti)q{t2)) = -- 9{ti -t2)+ dr ds g{ti - T)aR{T - s)g{t2 - s) 
I Jo Jo 

+ 9{ti)g{t2) {pD + g{h)g{h) (g^) + ( g{h)g{t2) + g{ti)g{h) ) ^ {q^p, + p.g.) , (48) 
c^R{r) = r duo {2n,{uj) + 1) . (49) 

The first term of the correlator comes from the anti-symmetric part, hence the com- 
mutator part [q{ti) , 5(^2) ], while the rest from the symmetric part. In this compu- 
tation the following identity was used; 



^2) = g{ti)g{t2) - g{ti)g{t2) +i du r{u) (h*{u, ti)h{u , ts) e-'^^*^"*^) - (c.c. 

Ju>c ^ 

(50) 
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Coincident time limits are evaluated from these, resulting in 

{q\t))= / du{2n,{u) + l)r{uj)\h{u,t)\^ 



+ 9\t) {VD + 9\t) {q-) + 9it)m {p,q, + q,p,) , (51) 
{p\t))= / du;{2ni{u) + l)r{u)\k{LJ,t)\^ 

+ g\t) {pI)+ 9^ it){q^) + g{t) 9 {t){p,q, + q,p,) , (52) 
- {q{t)p{t) + p{t)q{t)) = / dio {2ni{uj) + 1) r{uj) h{uj , t)k*{uj , t) 

+ 9{t)9{t) {pD + m '9 (t) (g5 + ( g\t) + g{t) 'g (t))^ {p,q, + q,p,) . (53) 
Asymptotic values are then 

(g2(oo)) = / dco ' ^ ' ^ ' {2n,ico) + 1) , (54) 



2uj 

{p\oo)) = duo- I (0|vl/(c.)) |2 {2n,{uj) + 1) . (55) 
More generally, for both t\ and ^2 in the asymptotic late time region, 

(g(ii)g(i2)) ^ / rfc^ ' ' cosL^(ti - ta) (2ni(L^) + 1) , (56) 



2cj 

|2 



kp{tr)v{t2)) ^ du-l {0\^itu)) \' cos^(ti - t2) {2ni{uj) + 1) . (57) 

The other correlator vanishes: {q{ti)p{t2) + p(t2)q(ti)) . 

As an illustration and for comparison with the result in the path integral, we 
take the initial environment in thermal state of temperature T = 1/(3 written by 
Q{uj) (and not diagonal Qiu)): 

(2n,(^) + l)=coth(^). (58) 

If we take for the initial system a gaussian state such as a thermal one having a 
different temperature from the environment (which includes the ground state in the 
limit of T — > 0), our results here coincide with previous formulas p|. 

In particular, we find it important to stress again that the asymptotic time limit 
of the occupation number, 

Moo) = -{^^^ + uq\cx^))-- 

= - du coth(^) {u + -) H{u) - - , (59) 
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with UJ the real part of the pole position, contains the fraction not suppressed by 
the usual Boltzmann factor. When the pole term dominates, or equivalently one 
approximates HiuJ^ by the Breit-Wigner function, then the temperature dependent 
part of the occupation number defined by 

n^=\ duj {uj + -) E{yj) , (60) 

has the factor e~'^/^ at low temperatures. But this approximation is not good at 
low temperatures. Indeed, let us examine a typical example by taking the form of 
r(co') = c{uj — cjc)" , with < a < 1 in the range oi uj^. < uJ <^ (^~^ cUc) and with 
UJ ^ Max ( ciJc , T ) . The result is 

~ iLr(a + l)e-^"=T°+^ (61) 

UJ'^ 

r is the Euler's gamma function. 

This shows that instead of the exponential suppression at low temperatures what 
is left in medium after the decay has a power-law behavior of temperature dependence 
(oc T""''^). An implication of this behavior to the unstable particle decay, as discussed 
in our previous paper 0, is that the remnant fraction in thermal medium does not 
suffer from the Boltzmann suppression factor at temperatures even much lower than 
the mass of the unstable particle. This is because the approximate Boltzmann-like 
equation is based on S-matrix elements computed on the mass shell, while the true 
quantum mechanical equation may contain quantities off the mass shell. Although 
implicit in these formulas (since tUc ~* limit is singular), the presence of a gap, 
ciJc = 2x (daughter mass) is critical to obtain the enhanced remnant fraction. 

The corresponding behavior of git) in this case is 

2c , , cos( cUct + 5 a ) , , 

,(t)^--r(a + l) ^ r (62) 

One can estimate the transient time from the exponential period to the power 
period by equating the two formulas of gii) in their respective ranges, to obtain 

( 2cr(afl)7-^0 ' ^''^ 

with 7 = — '^zq the decay rate. 

We finally mention briefly how the results derived by the operator solution are 
related to the path integral results. The gaussian system such as the one mentioned 



12 



is characterized by a gaussian reduced density matrix and is determined by the three 
quantities; (g^(t)) , , \ {q{t)p(t) + p{t)q{t)) . The results for these quantities 

given here agree with the results obtained by the path integral approach in our previ- 
ous work This shows that the two approaches are equivalent whenever both can 
give answers. We shall give more details on this equivalence in our forthcoming pa- 
per. The two approaches are complementary: the path integral approach can cope 
with more general dynamical systems than the harmonic one, while the operator 
approach, although limited to the linear harmonic system, offers a simple interpreta- 
tion in terms of the overlap probability of the initial system variable, with extention 
made possible to the general initial state. Both should be useful. 
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Figure Caption 



Figure 1 

Deformation of the complex contour for calculation of g{t). 
Figure 2 

Behavior of the function g(t) for the Ohmic spectrum weight r{uj) = c{uj — Uc) 
for Lo^ < u < Q. Chosen parameters are Uc = I , ^ = 100 , c = 0.01 , coq — 
2.11 . Only the positive region of g{t) > is shown for a technical reason, with clear 
signals of oscillation observed here. 
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